
COMPILE_OPT IDL2, LOGICAL_PREDICATE

conf = omni_load_conf()
s = omni_read_cat(conf.cat,ncat)
mw = omni_read_conffile('./conffiles/galactic_params.conf')
IF ~exist(pvec) THEN restore,'./local/'+conf.survey+'_pvec.sav',/ver
restore,'./local/'+conf.survey+'_velocities.sav',/ver

ind = where(constrain.kdist OR constrain.grsmatch OR $
            constrain.parallax OR constrain.knownd, ni)
p = pvec[ind]
c = constrain[ind]
s = s[ind]
v = v[ind]

print,ni



openw,lun,'./grs_paper/distcat.tex',/get_lun

ind = where(s.cnum GE 5500 AND s.cnum LE 5510, nl)

;; Select a good mix of sources to demonstrate the various properties
;;   of the table
ikn = where(p.stat.kdar EQ 'N', nkn)
ikf = where(p.stat.kdar EQ 'F', nkf)
ikt = where(p.stat.kdar EQ 'T', nkt)
iko = where(p.stat.kdar EQ 'O', nko)
iku = where(p.stat.kdar EQ 'U', nku)

nt = nkn+nkf+nkt+nko+nku
print,round(float(nkn)/nt*10.)>1,round(float(nkf)/nt*10.)>1,$
      round(float(nkt)/nt*10.)>1,round(float(nko)/nt*10.)>1,$
      round(float(nku)/nt*10.)>1

i1 = where(p.stat.kdar EQ 'N' AND v.rv_types EQ 16, n1)
i2 = where(p.stat.kdar EQ 'N' AND c.grsmatch, n2)
i3 = where(p.stat.kdar EQ 'F' AND v.rv_types EQ 1, n3)
i4 = where(p.stat.kdar EQ 'F' AND v.rv_types EQ 8, n4)
i5 = where(p.stat.kdar EQ 'T' AND c.grsmatch, n5)
i6 = where(p.stat.kdar EQ 'T' AND v.rv_types EQ 9, n6)
i7 = where(p.stat.kdar EQ 'O' AND v.rv_types EQ 9, n7)
i8 = where(p.stat.kdar EQ 'O' AND v.rv_types EQ 1, n8)
i9 = where(p.stat.kdar EQ 'U' AND c.grsmatch, n9)
i0 = where(p.stat.kdar EQ 'U' AND v.rv_types EQ 5, n0)

print,'Counts: ',n1,n2,n3,n4,n5,n6,n7,n8,n9,n0

i1 = i1[permute(n1)]
i2 = i2[permute(n2)]
i3 = i3[permute(n3)]
i4 = i4[permute(n4)]
i5 = i5[permute(n5)]
i6 = i6[permute(n6)]
i7 = i7[permute(n7)]
i8 = i8[permute(n8)]
i9 = i9[permute(n9)]
i0 = i0[permute(n0)]

ind = [i1[0],i2[0],i3[0],i4[0],i5[0],i6[0],i7[0],i8[0],i9[0],i0[0]]

print,'New i2: ',p[i2[0]].cnum
print,'New i5: ',p[i5[0]].cnum

ind = ind[sort(ind)]

;; ind = where_array([3283,3475,4969,5197,5560,6196,6729,6964,7612,8210],$
ind = where_array([3475,4061,4969,5197,5560,6196,6711,6964,7612,8210],$
                  p.cnum,count)
print,'COUNT:: ',count
nl = count

p = p[ind]
c = c[ind]
s = s[ind]
v = v[ind]
p.stat.dml /= 1.d3
p.stat.dbar /= 1.d3

FOR ii=0,nl-1 DO BEGIN
   
   ;; V2.1 Catalog Information
   str = string(s[ii].cnum, s[ii].glon LT 100 ? '\phn' : '',$
                s[ii].glon, s[ii].glat ge 0 ? '\phs': '',$
                s[ii].glat, s[ii].flux,s[ii].e_flux,$
                format="(I4,' & ',A0,F0.3,' & ',A0,'$',F0.3,'$ & $',F0.2,'(',F0.2,')$ & & ')")
   
   ;; Velocity Information
   IF v[ii].vlsr GE -500 THEN BEGIN
      IF v[ii].vlsr GE 0. THEN str += '\phs'
      IF abs(v[ii].vlsr) LT 100. THEN str += '\phn'
      IF abs(v[ii].vlsr) LT 10. THEN str += '\phn'
      str += '$'+string(v[ii].vlsr,format="(F0.1)")+'$ & '
      vlsr = v[ii].vlsr
      first = 1b
      FOR jj=0,6 DO BEGIN
         thisflag = v[ii].rv_types AND 2^jj
         IF thisflag THEN BEGIN
            IF ~first THEN str+=','
            CASE thisflag OF
               1:  str += '1'
               2:  str += '1'
               4:  str += '2'
               8:  str += (v[ii].l LT 120) ? '3' : '4' 
               16: str += '5'
            ENDCASE
            first = 0b
         ENDIF
      ENDFOR
   ENDIF ELSE BEGIN
      IF v[ii].grs.vlsr GE 0. THEN str += '\phs'
      IF abs(v[ii].grs.vlsr) LT 100. THEN str += '\phn'
      IF abs(v[ii].grs.vlsr) LT 10. THEN str += '\phn'
      str += '$'+string(v[ii].grs.vlsr,format="(F0.1)")+'$ & 6 '
      vlsr = v[ii].grs.vlsr
   ENDELSE
   str += ' & & '
   
   ;; Heliocentric Distance
   str += p[ii].stat.kdar+' & '+string(p[ii].stat.pml,format="(F0.2)")+' & '
   
   ;; Extract for ease of the user
   dml  = p[ii].stat.dml
   duse = p[ii].stat.duse
   dbar = p[ii].stat.dbar
   phstr = dml[0] LT 10. ? '\phn' : ''
   
   CASE p[ii].stat.kdar OF
      'N': str += phstr+string(dml[0],dml[1],dml[2],$
                         format="('$',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$')") + $
                  ' & \nodata & & '
      'F': str += phstr+string(dml[0],dml[1],dml[2],$
                         format="('$',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$')") + $
                  ' & \nodata & & '
      'O': str += phstr+string(dml[0],dml[1],dml[2],$
                         format="('$',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$')") + $
                  ' & \nodata & & '
      'T': str += phstr+string(dml[0],dml[1],dml[2],$
                         format="('$',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$')") + $
                  ' & '+ string(dbar[0],dbar[1],$
                                format="('$',F0.2,'(',F0.2,')$')") + '& & '
      'U': str += '\nodata & \nodata & & '
      ELSE:
   ENDCASE
   
   ;; Galactocentric Positions
   IF p[ii].stat.kdar EQ 'U' THEN BEGIN
      
      ell = p[ii].glon * !dtor
      ;; U sources, compute RGAL from VLSR
      rgal = (mw.v0 * sin(ell) * mw.r0 / (vlsr + mw.v0*sin(ell))) / 1.d3
      
      ;; Assume error of +- 7 km/s
      rgale1 = (mw.v0 * sin(ell) * mw.r0 / (vlsr -7.d + mw.v0*sin(ell))) / 1.d3
      rgale2 = (mw.v0 * sin(ell) * mw.r0 / (vlsr +7.d + mw.v0*sin(ell))) / 1.d3
      rgmeb = rgal - min([rgale1,rgale2])
      rgpeb = max([rgale1,rgale2]) - rgal
      rgeb = mean([rgmeb,rgpeb])
      
      phstr = rgal LT 10. ? '\phn' : ''
      str += phstr+string(rgal,rgeb,rgeb,$
                    format="('$',F0.2,'^{+',F0.2,'}_{-',F0.2,'}$')")+' & \nodata'
   ENDIF ELSE BEGIN
      
      ;; Non-U sources, compute RGAL from OMNI_LBD2RZ
      dists = [duse[0],duse[0]-duse[1],duse[0]+duse[2]]
      omni_lbd2rz,p[ii].glon,p[ii].glat,dists,Rgal,Z
      
      rgal /= 1.d3
      rgmeb = rgal[0] - min(rgal)
      rgpeb = max(rgal) - rgal[0]
      zmeb  = z[0] - min(z)
      zpeb  = max(z) - z[0]
      zeb = mean([zmeb,zpeb])
      
      ;; Add the systematic \pm 5 \kms from {Juric:2008} in quadrature.
      zeb = sqrt(zeb*zeb + 5.d*5.d)
      
      z_ph  = z[0] GE 0. ? '\phs' : ''
      z_ph += abs(z[0]) GE 10. ? '\phd' : '\phn'
      z_ph += abs(z[0]) LT 100. ? '\phn' : ''
      
      z_str   = abs(z[0]) GE 10. ? +string(round(z[0]),format="(I0)") : $
                +string(z[0],format="(F0.1)")
      zeb_str = zeb GE 10. ? string(round(zeb),format="(I0)") : $
                string(zeb,format="(F0.1)")
      
      phstr = rgal[0] LT 10. ? '\phn' : ''
      str += phstr+string(rgal[0],rgmeb,rgpeb,z_ph,z_str,zeb_str,format=$
                    "('$',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$ & "+$
                    "',A0,'$',A0,'(',A0,')$')")
      
   ENDELSE
   
   
   
   
   printf,lun,str+' \\'
   
ENDFOR

close,lun
free_lun,lun

;;******************************************************************************
;;******************************************************************************
;;******************************************************************************
;;******************************************************************************
;; Create MRT of Table 5

s = omni_read_cat(conf.cat,ncat)
restore,'./local/'+conf.survey+'_velocities.sav',/ver
ind = where(constrain.kdist OR constrain.grsmatch OR $
            constrain.parallax OR constrain.knownd, ni)
p = pvec[ind]
c = constrain[ind]
s = s[ind]
v = v[ind]

print,'PMRT: ',ni
help,s
print,m4_stat(ind)

openw,lun,'./grs_paper/table5.pmrt',/get_lun

FOR ii=0,ni-1 DO BEGIN
   
   ;; V2.1 Catalog Information
   str = string(s[ii].cnum, s[ii].glon, s[ii].glat, s[ii].flux,s[ii].e_flux,$
                format="(I4,' & ',F0.3,' & ',F0.3,' & ',F0.2,' & ',F0.2,' & ')")
   
   ;; Velocity Information
   IF v[ii].vlsr GE -500 THEN BEGIN
      str += ''+string(v[ii].vlsr,format="(F0.1)")+' & '
      vlsr = v[ii].vlsr
      first = 1b
      FOR jj=0,6 DO BEGIN
         thisflag = v[ii].rv_types AND 2^jj
         IF thisflag THEN BEGIN
            IF ~first THEN str+=','
            CASE thisflag OF
               1:  str += '1'
               2:  str += '1'
               4:  str += '2'
               8:  str += (v[ii].l LT 120) ? '3' : '4' 
               16: str += '5'
            ENDCASE
            first = 0b
         ENDIF
      ENDFOR
   ENDIF ELSE BEGIN
      str += ''+string(v[ii].grs.vlsr,format="(F0.1)")+' & 6 '
      vlsr = v[ii].grs.vlsr
   ENDELSE
   IF strmatch(str,'*1,1*') THEN BEGIN
      pos = strpos(str,'1,1')
      substr = strmid(str,0,pos) + strmid(str,pos+2)
      str = substr
   ENDIF
   
   str += ' & '
   
   ;; Heliocentric Distance
   str += p[ii].stat.kdar+' & '+string(p[ii].stat.pml,format="(F0.2)")+' & '
   
   ;; Extract for ease of the user
   dml  = p[ii].stat.dml / 1.d3
   duse = p[ii].stat.duse / 1.d3
   dbar = p[ii].stat.dbar / 1.d3

   CASE p[ii].stat.kdar OF
      'N': str += string(dml[0],dml[1],dml[2],$
                         format="(F0.2,' & ',F0.2,' & ',F0.2)") + $
                  ' & --- & --- & '
      'F': str += string(dml[0],dml[1],dml[2],$
                         format="(F0.2,' & ',F0.2,' & ',F0.2)") + $
                  ' & --- & --- & '
      'O': str += string(dml[0],dml[1],dml[2],$
                         format="(F0.2,' & ',F0.2,' & ',F0.2)") + $
                  ' & --- & --- & '
      'T': str += string(dml[0],dml[1],dml[2],$
                         format="(F0.2,' & ',F0.2,' & ',F0.2)") + $
                  ' & '+ string(dbar[0],dbar[1],$
                                format="(F0.2,' & ',F0.2)") + ' & '
      'U': str += '--- & --- & --- & --- & --- & '
      ELSE:
   ENDCASE
   
   ;; Galactocentric Positions
   IF p[ii].stat.kdar EQ 'U' THEN BEGIN
      
      ell = p[ii].glon * !dtor
      ;; U sources, compute RGAL from VLSR
      rgal = (mw.v0 * sin(ell) * mw.r0 / (vlsr + mw.v0*sin(ell))) / 1.d3
      
      ;; Assume error of +- 7 km/s
      rgale1 = (mw.v0 * sin(ell) * mw.r0 / (vlsr -7.d + mw.v0*sin(ell))) / 1.d3
      rgale2 = (mw.v0 * sin(ell) * mw.r0 / (vlsr +7.d + mw.v0*sin(ell))) / 1.d3
      rgmeb = rgal - min([rgale1,rgale2])
      rgpeb = max([rgale1,rgale2]) - rgal
      rgeb = mean([rgmeb,rgpeb])
      
      str += string(rgal,rgeb,rgeb,$
                    format="(F0.2,' & ',F0.2,' & ',F0.2)")+' & --- & ---'
   ENDIF ELSE BEGIN
      
      ;; Non-U sources, compute RGAL from OMNI_LBD2RZ
      dists = [duse[0],duse[0]-duse[1],duse[0]+duse[2]]
      omni_lbd2rz,p[ii].glon,p[ii].glat,dists*1.d3,Rgal,Z
      
      rgal /= 1.d3
      rgmeb = rgal[0] - min(rgal)
      rgpeb = max(rgal) - rgal[0]
      zmeb  = z[0] - min(z)
      zpeb  = max(z) - z[0]
      zeb = mean([zmeb,zpeb])
      
      ;; Add the systematic \pm 5 \kms from {Juric:2008} in quadrature.
      zeb = sqrt(zeb*zeb + 5.d*5.d)
      
      z_str   = abs(z[0]) GE 10. ? string(round(z[0]),format="(I0)") : $
                string(z[0],format="(F0.1)")
      zeb_str = zeb GE 10. ? string(round(zeb),format="(I0)") : $
                string(zeb,format="(F0.1)")
      
      str += string(rgal[0],rgmeb,rgpeb,z_str,zeb_str,format=$
                    "(F0.2,' & ',F0.2,' & ',F0.2,' & "+$
                    "',A0,' & ',A0)")
      
   ENDELSE
   
   
   
   
   printf,lun,str;+' \\'
   
ENDFOR

close,lun
free_lun,lun

;; Write out PMRT header info for use with AAS tool thingie
openw,lun,'./grs_paper/table5.hdr',/get_lun

printf,lun,'---'
printf,lun,'deg'
printf,lun,'deg'
printf,lun,'Jy'
printf,lun,'Jy'
printf,lun,'km/s'
printf,lun,'---'
printf,lun,'---'
printf,lun,'---'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'kpc'
printf,lun,'pc'
printf,lun,'pc'

printf,lun,''
printf,lun,''
printf,lun,''

printf,lun,'CNUM'
printf,lun,'GLON'
printf,lun,'GLAT'
printf,lun,'Flux'
printf,lun,'e_Flux'
printf,lun,'VLSR'
printf,lun,'r_VLSR'
printf,lun,'KDAR'
printf,lun,'PML'
printf,lun,'dML'
printf,lun,'e_dML'
printf,lun,'E_dML'
printf,lun,'dbar'
printf,lun,'e_dbar'
printf,lun,'Rgal'
printf,lun,'e_Rgal'
printf,lun,'E_Rgal'
printf,lun,'z'
printf,lun,'e_z'

printf,lun,''
printf,lun,''
printf,lun,''

printf,lun,'Bolocat V2.1 Catalog Number'
printf,lun,'Galactic Longitude'
printf,lun,'Galactic Latitude'
printf,lun,'Integrated BGPS Flux Density'
printf,lun,'Standard Error on Integrated Flux Density'
printf,lun,'LSR Velocity'
printf,lun,'Velocity Reference (1)'
printf,lun,'KDA Resolution (2)'
printf,lun,'Maximum-likelihood Probability (3)'
printf,lun,'?Maximum-likelihood distance'
printf,lun,'?Lower error in dML'
printf,lun,'?Upper error in dML'
printf,lun,'?Weighted-average distance'
printf,lun,'?Weighted-average distance error'
printf,lun,'Galactocentric Radius'
printf,lun,'Lower error in Rgal'
printf,lun,'Upper error in Rgal'
printf,lun,'?Galactocentric Vertical Position (4)'
printf,lun,'?Galactocentric Vertical Position error (5)'

close,lun
free_lun,lun

END
